Lecture (Week 11b): Clustering
1. Data
data(faithful)
dim(faithful)[1] 272 2with(faithful, plot(eruptions, waiting, main="Old Faithful Geyser"))
2. \(K\)-means
(r0 = kmeans(faithful, centers=2)) # K = 2K-means clustering with 2 clusters of sizes 172, 100
Cluster means:
eruptions waiting
1 4.29793 80.28488
2 2.09433 54.75000
Clustering vector:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
1 2 1 2 1 2 1 1 2 1 2 1 1 2 1 2 2 1 2 1
21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
2 2 1 1 1 1 2 1 1 1 1 1 2 1 1 2 2 1 2 1
41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
1 2 1 2 1 1 2 2 1 2 1 1 2 1 2 1 1 2 1 1
61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80
2 1 2 1 2 1 1 1 2 1 1 2 1 1 2 1 2 1 1 1
81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
1 1 1 2 1 1 1 1 2 1 2 1 2 1 2 1 1 1 2 1
101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
2 1 2 1 1 2 1 2 1 1 1 2 1 1 2 1 2 1 2 1
121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140
2 1 1 2 1 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1
141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160
1 2 1 1 1 2 1 2 1 2 1 1 2 1 1 1 1 1 2 1
161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180
2 1 2 1 2 1 2 1 2 1 2 2 1 1 1 1 1 2 1 1
181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200
2 1 1 1 2 1 1 2 1 2 1 2 1 1 1 1 1 1 2 1
201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220
2 1 1 2 1 2 1 1 2 1 1 1 2 1 2 1 2 1 2 1
221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240
2 1 2 1 1 1 1 1 1 1 1 2 1 2 1 2 2 1 1 2
241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260
1 2 1 2 1 1 2 1 2 1 2 1 1 1 1 1 1 1 2 1
261 262 263 264 265 266 267 268 269 270 271 272
1 1 2 1 2 2 1 1 2 1 2 1
Within cluster sum of squares by cluster:
[1] 5445.591 3456.178
(between_SS / total_SS = 82.4 %)
Available components:
[1] "cluster" "centers" "totss" "withinss" "tot.withinss"
[6] "betweenss" "size" "iter" "ifault" with(faithful, plot(eruptions, waiting, col=r0$cluster+1, main="K = 2 (Original Data)"))
points(r0$centers, col=2:3, pch=19, cex=2)
Should we standardise the data first?
faithful2 = as.data.frame(scale(faithful)) # then each variable has mean 0 and standard deviation 1
head(faithful2) eruptions waiting
1 0.09831763 0.5960248
2 -1.47873278 -1.2428901
3 -0.13561152 0.2282418
4 -1.05555759 -0.6544374
5 0.91575542 1.0373644
6 -0.52987412 -1.1693335attach(faithful2)
r = kmeans(faithful2, centers=2) # K = 2
par(pty="s") # square plotting region
plot(eruptions, waiting, col=r$cluster+1, main="K = 2 (Standardised Data)")
points(r$centers, col=2:3, pch=19, cex=2)
cex = 1.5
par(mfrow=c(2,2), pty="s")
for(k in 3:6) {
r = kmeans(faithful2, centers=k)
plot(eruptions, waiting, col=r$cluster+1, main=paste0("K = ", k),
cex=cex, cex.axis=cex, cex.lab=cex, cex.main=cex)
points(r$centers, col=2:(k+1), pch=19, cex=2)
}
3. Adjusted Rand Index:
library(mclust)Package 'mclust' version 5.4.7
Type 'citation("mclust")' for citing this R package in publications.ari = double(5)
for(k in 2:6) {
r = kmeans(faithful2, centers=k)
ari[k-1] = adjustedRandIndex(r0$cluster, r$cluster)
}
ari[1] 0.9415368 0.5789191 0.4556733 0.3459902 0.30580344. Mixture-based Clustering
(r = Mclust(faithful, G=2, modelNames="VVV"))'Mclust' model object: (VVV,2)
Available components:
[1] "call" "data" "modelName" "n"
[5] "d" "G" "BIC" "loglik"
[9] "df" "bic" "icl" "hypvol"
[13] "parameters" "z" "classification" "uncertainty" summary(r)----------------------------------------------------
Gaussian finite mixture model fitted by EM algorithm
----------------------------------------------------
Mclust VVV (ellipsoidal, varying volume, shape, and orientation) model with 2
components:
log-likelihood n df BIC ICL
-1130.264 272 11 -2322.192 -2322.697
Clustering table:
1 2
175 97 plot(r, "uncertainty")
title(main="K = 2")
predict(r)$classification # clustering [1] 1 2 1 2 1 2 1 1 2 1 2 1 1 2 1 2 2 1 2 1 2 2 1 1 1 1 2 1 1 1 1 1 1 1 1 2 2
[38] 1 2 1 1 2 1 2 1 1 1 2 1 2 1 1 2 1 2 1 1 2 1 1 2 1 2 1 2 1 1 1 2 1 1 2 1 1
[75] 2 1 2 1 1 1 1 1 1 2 1 1 1 1 2 1 2 1 2 1 2 1 1 1 2 1 2 1 2 1 1 2 1 2 1 1 1
[112] 2 1 1 2 1 2 1 2 1 2 1 1 2 1 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 1 2 1 1 1 2 1 2
[149] 1 2 1 1 2 1 1 1 1 1 2 1 2 1 2 1 1 1 2 1 2 1 2 2 1 1 1 1 1 2 1 1 2 1 1 1 2
[186] 1 1 2 1 2 1 2 1 1 1 1 1 1 2 1 2 1 1 2 1 2 1 1 2 1 2 1 2 1 1 1 2 1 2 1 2 1
[223] 2 1 1 1 1 1 1 1 1 2 1 2 1 2 2 1 1 2 1 2 1 2 1 1 2 1 2 1 2 1 1 1 1 1 1 1 2
[260] 1 1 1 2 1 2 2 1 1 2 1 2 1par(mfrow=c(2,2))
for(k in 3:6) {
r = Mclust(faithful, G=k, modelNames="VVV")
plot(r, "uncertainty",
cex=cex, cex.axis=cex, cex.lab=cex)
title(main=paste0("K = ", k), cex.main=cex)
}
5. Hierarchical Clustering
Let’s still use standardised data.
d = dist(faithful2) # pairwise Euclidean distances
class(d)[1] "dist"head(d)[1] 2.4225392 0.4358752 1.7014945 0.9289700 1.8737969 1.1692204Complete linkage
r = hclust(d) # complete linkage, by default
names(r)[1] "merge" "height" "order" "labels" "method"
[6] "call" "dist.method"plot(r, cex.axis=cex, cex.lab=cex, cex.main=cex) # dendrogram
cutree(r, 2) # cluster labels of observations for K = 2 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
1 2 1 2 1 2 1 1 2 1 2 1 1 2 1 2 2 1 2 1
21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
2 2 1 1 1 1 2 1 1 1 1 1 1 1 1 2 2 1 2 1
41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
1 2 1 2 1 1 1 2 1 2 1 1 2 1 2 1 1 2 1 1
61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80
2 1 2 1 2 1 1 1 2 1 1 2 1 1 2 1 2 1 1 1
81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
1 1 1 2 1 1 1 1 2 1 2 1 2 1 2 1 1 1 2 1
101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
2 1 2 1 1 2 1 2 1 1 1 2 1 1 2 1 2 1 2 1
121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140
2 1 1 2 1 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1
141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160
1 2 1 1 1 2 1 2 1 2 1 1 2 1 1 1 1 1 2 1
161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180
2 1 2 1 1 1 2 1 2 1 2 2 1 1 1 1 1 2 1 1
181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200
2 1 1 1 2 1 1 2 1 2 1 2 1 1 1 1 1 1 2 1
201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220
2 1 1 2 1 2 1 1 2 1 2 1 2 1 1 1 2 1 2 1
221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240
2 1 2 1 1 1 1 1 1 1 1 2 1 2 1 2 2 1 1 2
241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260
1 2 1 2 1 1 2 1 2 1 2 1 1 1 1 1 1 1 2 1
261 262 263 264 265 266 267 268 269 270 271 272
1 1 2 1 2 2 1 1 2 1 2 1 par(mfrow=c(2,2))
for(k in 2:5) # complete linkage
plot(eruptions, waiting, col=cutree(r,k)+1, main=paste0("K = ", k),
cex=cex, cex.axis=cex, cex.lab=cex, cex.main=cex)
Single linkage
par(mfrow=c(2,2))
r = hclust(d, method="single") # single linkage
for(k in 2:5) # single linkage
plot(eruptions, waiting, col=cutree(r,k)+1, main=paste0("K = ", k),
cex=cex, cex.axis=cex, cex.lab=cex, cex.main=cex)
Average linkage
par(mfrow=c(2,2))
r = hclust(d, method="average") # average linkage
for(k in 2:5) # average linkage
plot(eruptions, waiting, col=cutree(r,k)+1, main=paste0("K = ", k),
cex=cex, cex.axis=cex, cex.lab=cex, cex.main=cex)
Centroid linkage
par(mfrow=c(2,2))
r = hclust(d, method="centroid") # centroid linkage
for(k in 2:5) # average linkage
plot(eruptions, waiting, col=cutree(r,k)+1, main=paste0("K = ", k),
cex=cex, cex.axis=cex, cex.lab=cex, cex.main=cex)
detach(faithful2)6. Heatmap
heatmap(as.matrix(faithful), scale="column", distfun=dist, hclustfun=hclust,
margins=c(15,5))